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ABSTRACT 

Our understanding of the mechanisms governing the structure and secular evolution 
galaxies assume nearly integrable Hamiltonians with regular orbits; our perturbation 
theories are founded on the averaging theorem for isolated resonances. On the other 
hand, it is well-known that dynamical systems with many degrees of freedom are ir¬ 
regular in all but special cases. For example, the solar system is chaotic and no doubt 
more so when dynamically young. Similarly, although driven towards symmetry by 
dissipation and mixing galaxies are chaotic and, galaxies are dynamically young, are 
strongly affected by perturbations, and are probably strongly chaotic. The best de¬ 
veloped framework for studying the breakdown of regularity and the onset is the 
Kolmogorov-Arnold-Moser (KAM) theory. Here, we use a numerical version of the 
KAM procedure to construct regular orbits (tori) and locate irregular orbits (bro¬ 
ken tori). Irregular orbits are most often classified in astronomical dynamics by their 
exponential divergence using Lyapunov exponents. Although their computation is nu¬ 
merically challenging, the procedure is straightforward and they are often used to esti¬ 
mate the measure of regularity. The numerical KAM approach has several advantages 
over Lyapunov exponent computation: 1) it provides the morphology of perturbed 
orbits; 2) its constructive nature allows the tori to be used as basis for studying sec¬ 
ular evolution; 3) for broken tori, clues to the cause of the irregularity may be found 
by studying the largest, diverging Fourier terms; and 4) it is more likely to detect 
weak chaos and orbits close to bifurcation. Conversely, it is not a general technique 
and works most cleanly for small perturbations. We develop a perturbation theory 
that includes chaos by retaining an arbitrary number of interacting terms rather than 
eliminating all but one using the averaging theorem. The companion papers show that 
models with significant stochasticity seem to be the rule, not the exception. Moreover, 
the subdimensionally of phase-space features spawned by chaotic dynamics demands 
fine-scale temporal resolution, and we discuss the implication for n-body simulations. 

Key words: methods: numerical — dark matter — galaxies: formation, evolution, 
kinematics and dynamics — chaos 


1 INTRODUCTION 
1.1 Overview 

Much of our insight into the dynamics of galaxy structure and evolution relies on the properties of regular orbits. We commonly 
assume that regular orbits dominate owing to the importance of baryonic dissipation during galaxy formation: dissipation 
naturally selects non-intersecting orbits and is the principal enforcer of axisymmetry. Following the mostly successful ACDM 
paradigm, we assume that galaxies are a combination of dissipative (disks) and non-dissipative (dark-matter and stars) 
components with different geometries. There is no natural reason to assume that sufficient conditions for regularity (i.e. the 
simultaneous separability of both the Hamilton-Jacobi equation and the Poisson equation defined by the Stackel conditions. 
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e.g. Vinti, Der & Bonavito 1998) coincides with this combination of dissipative and dissipationless dynamical processes. 
Therefore, these dynamical systems are not regnlar everywhere almost certainly. 

On the other hand, owing to ease of computation, insight provided by regular systems underpins most of our physical 
scenarios for galactic stability, structure and evolution. For example, theories of spiral arms, swing amplification, among 
others assume regularity. The author has made predictions for the secular evolution of galaxies in combined disk and spheroid 
systems using classical Hamiltonian perturbation theory (e.g. Tremaine & Weinberg 1984; Weinberg 1998, 2001; Vesperini 
& Weinberg 2001; Weinberg & Katz 2002). This previous work assumes that only one resonance influences the dynamics of 
any one otherwise regular orbit at one time using the averaging theorem (e.g. Arnold 1989). The current paper assesses the 
failure of this assumption and provides a tool for identifying important dynamics missed by the averaging approximation. 
Moreover, environmental perturbations and secular and possibly dynamical instabilities yield bars and spiral arms. Some of 
these features are strong perturbations and are likely to generate irregular orbits. Specifically, is it possible that our ignorance 
of irregular behaviour biases our interpretive insight? 

In this paper, we present a general method based on one of the mathematical tools used in some proofs of the Kolmogorov- 
Arnold-Moser (KAM) theorem for studying the regular orbits that are preserved under perturbations. Recall that regular 
orbits in a bound system are fully characterised by a full set of constants of the motion, often chosen to be action-angle variables 
for theoretical convenience. These variables are constructed by finding the canonical transformation from a set of Cartesian 
variables (or some other conical coordinate system) to a new set which leaves the Hamiltonian independent of angles. The 
Hamilton-Jacobi (HJ) equation defines this canonical transformation implicitly. For regular systems, the action-angle variables 
define tori in phase space that fully describe the quasi-periodic motion. The approach pursued here is the construction of 
numerical solution of the HJ equation for a perturbed regular system. If the solution exists, the new transformation yields a 
new set of action-angle variables. If the solution does not exist, the torus is presumed destroyed; that is, one or more of the 
actions is no longer a constant of the motion and the trajectory may appear chaotic. 

This approach is complementary to the most popular approach for studying chaos: Lyapunov exponent analysis. While a 
study of Lyapunov exponents provides a local phase-space measure for orbits whose trajectories diverge exponentially under 
the influence of a perturbation, exponent values close to zero provide little useful information for the surviving regular orbits. 
Conversely, the KAM procedure characterises the orbits distorted by a perturber while identifying but not characterising the 
orbit of the broken tori. Finally, if irregular dynamics is important for galaxies, the dynamics are likely to appear in stochastic 
layers around primary resonances and these need to be accounted for by criteria used by n-body simulators to determine 
phase-space and temporal resolution (see Springel 2005 for a widely-used example). This paper develops the basic formalism 
and presents the results for the circular restricted three-body problem as an example. In the companion papers, we explore 
the question of how prevalent is irregularity in systems of interest and how might this irregularity affect galaxy structure and 
evolution? 

Since these issues in modern non-linear dynamics will be unfamiliar to many, we begin with a quick review of the 
theoretical motivation in section 1.2 followed by a discussion of standard practice in galaxy dynamics section 1.3. 


1.2 Underlying dynamical theory 

KAM theory is at the heart of modern Hamiltonian dynamics. The fundamental problem addressed by KAM theory first 
arose in the three-body problem beginning with Poincare: the introduction of a third body into the integrable gravitational 
two-body system renders the dynamics non-integrable, independent of the mass of the bodies. There exists a large literature 
on perturbative tools for constructing approximate integrable systems for solar-system dynamics (see Brouwer & Clemence 
1961, for a pre-KAM summary). These tools are based on expansions the ratio of the planet mass to the solar mass, a small 
parameter. Then, they reduce the degrees of freedom of the problem by identifying uninteresting degrees of freedom and 
integrating or averaging the resulting equations of motion over these uninteresting variables. However, the convergence of 
these series is questionable owing to the appearance of small denominators in the solution to the averaged equations. Poincare 
called this fundamental problem of mechanics (see Barrow-Green 1997 for an historical review). 

The essence of the vanishing-denominator conundrum will be familiar to most physics students in the context of a forced 
harmonic oscillator. If the external force matches the oscillator’s natural frequency (i.e. is resonant) the linear solution is 
unbounded. More generally, the response of the non-linear oscillator is bounded but will tend to generate oscillations in 
multiples of the driving frequency resulting many resonances. If we add a second non-linear oscillator, the accidental near 
coincidence of two harmonics may lead to erratic behaviour. In a similar way, let us consider a solar system with two planets; 
one planet exerts a periodic force on the motion of a second. Expanding the perturbation in a Fourier series, we may recast 
the problem as that of coupled pendula forced by all combinations of all possible harmonics of the natural frequencies of the 
two planets. The commensurability of the orbital periods for one or more terms in the Fourier series of sufficient amplitude 
can lead to small denominators resulting from resonant and divergence of the expansion. 

Kolmogorov (1954) suggested an two-part approach to obtaining convergence of these series. First, proceeding in the 
standard way, one may linearize the problem about an approximate solution and obtain an the approximate solution complete 
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with small denominators. Secondly, one may improve the approximate solution inductively using a qnadratically convergent 
iterative method. Success of this approach relies on the choice of initial conditions such that the denominators do not vanish 
for low-order harmonic terms and that the amplitude of the harmonics terms approach zero exponentially with increasing 
harmonic order. These ideas were then made rigorous by Arnold (i.e. 1963) and Moser (1962, 1966) and led to what we now 
know as the KAM theory. In short, systems that can be shown to satisfy the KAM theorem have to be snfficient smooth and 
continuous and not have commensurate frequencies to some maximal order for some perturbation amplitude e smaller than 
some bonnd, e < e*. 

KAM theory does not solve the problem of non-integrability but helps explain its genesis. Specifically, tori may be 
destroyed if the frequencies approach low-order commensurabilities or the perturbation strength e grows. Although KAM 
theory does not estimate the chaotic measure, direct numerical solution of the perturbed equations of motion or evolution 
using Hamiltonian maps suggests that modest perturbation strengths induce a non-zero measure of chaotic orbits. The 
resulting phase space consists of a combination of regular and irregular regions often called a mixed phase space. When most 
tori are broken, even those far from the original location of the vanishing denominators, the actions of the chaotic trajectories 
drift in the irregular regions and the remaining tori are localised to resonances, so-called regular islands. A detailed description 
of this theory is beyond the scope of this paper. However, while reading various introductions to the current version of the 
KAM theory (e.g. De La Llave 2003) while on sabbatical at the Institute for Advanced Study several years ago, I was struck 
by the algorithmic nature of Kolmogorov’s underlying concept. The basic tools are part of the every numericist’s standard 
arsenal: Newton’s method and numerical Fourier transforms. Can one use these tools numerically to understand practical 
astronomical problems? The answer appears to be: YES! The outstanding question for us here is: does this generic dynamical 
picture of mixed phase space have implications for galactic dynamics? 


1.3 Application to galactic dynamics 

The same general ingredients necessary to understand modern celestial dynamics applies to our understanding of galactic 
dynamics. As in the perturbation theories described in Tremaine & Weinberg (1984); Weinberg (1994), we may exploit 
the quasi-periodic nature of the unperturbed representation of the phase space to convert the dynamical partial-differential 
equations to an series of algebraic equations. This series corresponds to an ordering in increasing spatial frequencies. Truncating 
this series on physical grounds yields an algebraically tractable set of linear equations. The main difference arises in the broader 
frequency spectrum and therefore slower convergence the perturbed Kepler problem. In addition, the approach here abandons 
the averaging theorem that artihcially isolated resonances. 

Using a related but orthogonal approach, James Binney and collaborators have made enormous progress in their program 
to derive approximate angle-action coordinates from the standard Cartesian phase-space coordinates (and vice versa) that 
they call torus mapping. In a fiducial series of papers, McGill & Binney (1990); Binney & Kumar (1993); Kaasalainen & Binney 
(1994) show that generating functions formulated as action-angle series together with a careful choice of a nearly analytic 
Hamiltonian that approximates the system of interest and an appropriately tuned cost function may be used to solve for 
the coefficients of the canonical generating function with sufficient accuracy to obtain tori that are excellent approximations 
to systems of astronomical interest. More recently, Sanders & Binney (2014) demonstrated that this approach works well in 
three-dimensional as well as two-dimensional problems. 

Although both approaches begin with an action-angle expansion of the canonical generating function, we are attempting 
to discover, quantify and understand the role of irregularity galaxies rather than trying to obtain regular approximations to 
galaxies. Torus mapping uses a cost function to match the coefficients of the generating function series in the least-squares 
sense to a target Hamiltonian. In the approach considered here, we attempt to implicitly solve for the values of these coefficients 
that leaves the target Hamiltonian invariant. Formally, this is equivalent to solving the Hamilton-Jacobi equation that defines 
the actions and angles, as described earlier. 

As a first galactic application, Weinberg (2015a, hereafter Paper 2) applies the method developed in this paper to 
investigate the role of chaos in barred galaxies. Most certainly, we are not the first to consider the importance of chaos in 
barred galaxies; e.g. for a selection of relatively recent papers see Patsis (e.g. 2012); Manos, Bountis & Skokos (e.g. 2013); 
Contopoulos & Harsoula (e.g. 2013); Manos & Machado (e.g. 2014a); Ernst & Peters (e.g. 2014) and references therein. 
Rather, this work was intended to test the approach on a well investigated astronomical scenario. Using the new numerical 
KAM approach, we can easily identify, in phase space, the tori that survive under the self-consistent perturbation of the bar 
potential on the overall gravitational potential. 

Secondly, although collisionless and collisional dissipation lead to systems that are approximately regular, the combination 
of the two are not regular. Specifically, a dark halo yields a weakly triaxial system, and therefore approximately spherical 
system. Stable spherical systems are regular. The baryonic component quickly settles into an approximately axisymmetric or 
elliptical disk, another nearly regular system. However the join of the two systems, as we envision most galaxies, are unlikely 
to be regular. Moreover, a general three-dimensional stellar disk will not have regular orbits. Therefore, there will be irregular, 
chaotic regions in phase space in addition to regular regions. This simple observation raises several questions. Are there generic 
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features of a exponential stellar disk induced by breaking of invariant tori? Do the chaotic regions enable significant structural 
evolution in a galaxy lifetime? Therefore, as a second application of the numerical KAM method, we consider this problem in 
Weinberg (2015b, hereafter Paper 3). 


1.4 Plan for this paper 

To investigate these questions, we will develop a computational method that seeks numerical solutions to the Hamilton- 
Jacobi (HJ) equation for a irregular system by successive perturbations to an initially regular system. We will review the 
mathematical background and develop the method in section 2. For the galactic applications described above, we will begin 
with the spherical system whose radial profile is obtained from the monopole or other appropriate spherical approximation 
for gravitational potential of the full non-spherical mass distribution. The perturbation, then, is the difference between full 
mass distribution and the monopole distribution. For example, to investigate irregularity of the combined dark halo-galactic 
disk system, we consider the unperturbed system to be the monopole part of the multipole expansion for the entire system. 
The resulting massless perturbation squeezes the isopotential surfaces towards the disk plane. Then, one may investigate the 
evolving regular orbits (tori) by successively increasing the perturbation until the target mass distribution is reached, section 
1.2 describes the mathematical approach followed by a description of the numerical algorithm whose details can be found in 
Appendix A. We conclude with a summary and discussion in section 6. 


2 USING THE KAM IDEAS NUMERICALLY 
2.1 Introduction 


The action-angle variables for a time-independent equilibrium system are the result of seeking a canonical transformation that 
leaves the Hamiltonian a function of momentum variables alone. If the ansatz of a generating function separable in the new 
momentum variables also separates the Hamiltonian for bound values of the energy, the resulting the resulting equations of 
motion are angle-like variables that increase linearly with time describing the quasi-periodic motion. The new momenta are 
constant in time because the Hamiltonian is coordinate independent. With an appropriate scaling, the new coordinates may 
be chosen to advance by 27r for each cycle. Let these new momenta and coordinates (action-angle variables) be denoted by 
(I, w). Geometrically, these new coordinates define a tori in phase space. If the system is everywhere regular, the phase space 
is foliated by tori. 

The existence of action-angle variables implies a system of regular orbits. Although this construction is sufficient but 
not necessary for regularity, it is the usual starting point for most discussions. More generally, orbits are often defined to 
be irregular if distance between two phase-space states with small initial distance increases exponentially. This exponential 
divergence implies that any notion of predictability or memory of the initial state is soon lost. Clearly, trajectories defined by 
action-angle variables are confined to tori and will not be exponentially divergent. 

Now, assume a system of regular orbits specified by action-angle variables (I, w) with the Hamilton Ho (I). For the 
problems considered in this paper, we will begin with a spherical unperturbed mass distribution. In the case, the values for 
(I, w) are always available by quadrature given Ho (I). Thus, the regular system defined by a spherical potential provides a 
good starting point for adding perturbations that induce irregularity. The same features hold true for any potential in Stackel 
form (Stackel 1892) but the general coordinate systems are more cumbersome (e.g. Binney & Tremaine 2008). As we will see 
in the development below, a chain of successive perturbations and canonical transformations describing invariant tori will be 
tied, through interpolation, to this fiducial foliation of phase space by regular orbits. We apply a non-spherical perturbation. 
Hi (I, w) = Vi(r, 9, <j>) where r, 9, 0 are the radius, colatitude, and azimuthal angle, respectively, so that the total Hamiltonian 
is 


H(I,w) = Ho(I)-fHi(I,w). 


( 1 ) 


Given some non-spherical or non-axisymmetric perturbation, we often rewrite equation (1) to define the perturbation: 
Hi(I,w) = H(I, w) — Ho(I). For a spherical distribution I and w are three-vectors (see Tremaine & Weinberg 1984, for 
an explicit representation). We could use the extended phase-space paradigm (e.g. Wiggins 2003, Chap. 3) to include time de¬ 
pendence without any substantive changes to this development. However, to simply notation, we assume that the perturbation 
is time independent. 

Following the standard procedure for solving the HJ equation, we seek a new set of canonical variables, (J,s) say, such 
that H„eu)(J) = H(I,w). We assume a time-independent canonical generator to transform from the original action-angle 


variables to a new set of action-angle variables (I, w) 


Fi (J,w) 


■> (J, s) to first order in the perturbation: 


Fi(J, w) = J • w -I- Wi{J, w). 


( 2 ) 


This is a standard Type 2 generating function (Goldstein, Poole & Safko 2002). The first term in equation (2) generates the 
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identity transformation. The new canonical variables are then (e.g. Arnold 1989, §48) 


and the new Hamiltonian is 


I = 

s = 


dFi 

9w 

dFi 

1)1 


= J + 


= w + 


dWi 

9w 

dWi 

HT 


[old] 

[new] 


J + 


dW^ 
dw ' 


Ho{J + 



+ Hi 



dWi 

9w 



Recall that 


(3) 

(4) 

(5) 


= ( 6 ) 

defines the time-independent HJ equation and is an implicit equation for the term Wi through equation (3). The key to 
finding is the implicit solution for Wi defined by equation (5). Once we have the generating function Wi, the new 

action variables follow immediately from equation (3). 


2.2 The iterative canonical perturbation theory 

A direct solution of the HJ equation (defined through equations 5 and 6) using a numerical partial differential equation scheme 
is not practical for problems of interest. Rather, we define a solution based on a succession of first-order expansions. The 
solution is then iterated to improve the approximation to Wi until the desired convergence is obtained. 


2.2.1 Derivation of the algorithm 

The solution is constructed inductively as follows. Assume that we have an approximate solution w["^ at iteration n whose 
deviation from an exact solution is 

q^"^ = H(^J+ - H^L (J) (7) 

where 

(J) = ^ J + ^, w) . (8) 

The solution will be OHi, that is, the order of the perturbing potential. If Wl"^ is exactly Wi, then equation (7) yields 
= 0 owing to equation (5). Similarly, if is exactly Hnew{J) then the definition in equation (8) is an identity. 

However, if ^ 0, then is second-order in the small perturbation by construction. Now, let us look for an improved 
solution, -|- where C1[A1T|"^] = Cl[(/f”^] from the definition in equation (7). The improvement AWl"^ 

removes an additional piece of the discrepancy between hIZI and Hnew which is already of order We will explicitly 
derive the improvement in terms of the next order difference gN+i] ^nd demonstrate the following ordering 

qln+l] =h(^J+ = 0[(q'"')"]. (9) 

This implies that the iteration converges quadratically. We assume the existence of solution and endeavour to obtain 
the solution such that 

(J) = Hi^J+ , wj . (10) 


Although we do not know (J) a priori, it is a constant with respect to the new angle s will not affect our solution. We 

begin with the Taylor series expansion of equation (10) in AW]"': 


HtllA) = h(j+ =h(3+ -f 


i9w 


dWC 

= H-[J+^^,w| + 

O'W 


dw 

dH dAW["^ 


dl 


9w 


+ 0[{q 


d'w 


,w 


= g'"’ + M"l(J) + ^ + o[(g'"’)1 


( 11 ) 
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having substituted equation (7) and (dH/dl) is evaluated at I = J + /d'w). The function {dH/d\) may be evaluated 

implicitly in terms of iteration quantities as follows: 


dq'"' 

9J 


dH 


d3 


9J 


dH 
dl ' 

1 + 


1 + 


9J9w 


After rearranging, we find 


Since 


dH 

IW 


1 + 


dwp 

9Ji9w 


1 + 


d^wp 

9J9w 


a?'" 

aj 


d^wp 

dJdw 


aj 


aAiTj” 

aw 


aj 

dH dHPpj) 

ai 


aj 


dnP^ 

dj 


= o 




( 12 ) 


(13) 


(14) 


we may ignore the dq^^^/dJ term in equation (13) upon substituting into equation (11) to 0[{q^"^p. Performing the substi¬ 
tution, the equality in the last of equation (11) becomes 

rM 


dH^ 

dJ 


1 + 


up to terms with 0[{q^"P\. Because (p\ hP-u, and are known quantities at this point in the iteration, equation (15) 
may be solved explicitly for AwP to first-order accuracy yielding a solution to = 0 that is accurate to Cl[(A(/f”^)^]. In 

other words, we have shown that = Cl[(Agf"^)^], and, therefore, the iterative correction to defined by equations 

(7) and (15) converges quadratically. 

A compact solution of equation (15) is obtained by Fourier analysis in the new angles s using the standard Type 2 
canonical transformation (equation 4). We assume that the action-angle expansion of the generating function Wi in the new 
action-angle variables exists: 

Wi = ^r!7i(J)e”''". (16) 


d^wp 

dSd'w 


d^WP _ 

dw 


rrl'H-il/'T'i _ /-TN 

-^-^new V / -^-^new\^J 


(15) 


The Fourier transform of equation (15) may be written 

qi = y^^;^^dsg'’"'(J,w(J,s))e 


(2^) 

1 


ds 


dHil 


dJ 


(27r)3 

-( HpPip-nPpj)^ 


1 + 


dpvp 

dJdw 

-il-s 


dAWp 

dw 


(17) 


The first term in the integrand of equation (17) may be simplified using the canonical transformation (equation 4): 


1 -f 


d- 




dJdw J 


OAlTi 

dw 


1-f 

9s 

dw 


d^Wi 


dwdJ 
1 


dAWi 


dw 

dAWi dAWi 


(18) 


9w ds 

The second term in the integrand is independent of s and therefore vanishes unless 1 = 0. Finally, in terms of the Fourier 
coefficients, the solution to equation (15) is 

dnP 


or, for 1/0, 


qi + {HPL - HpPpio + il ■ 


Arc7i(J) = 


9J 


-Arc7i = 0 


1- 


SHpl 

dJ 


-q\- 


(19) 


( 20 ) 


Since an arbitrary additive constant in the generating function W\ does not affect the solution for the action-angle variables, 
we are free to ignore the constant coefficients with 1 = 0. Repeating the procedure, we may continue to iterate until q^"^ is as 
small as desired. This iteration method is essentially that used in proving the KAM theorem (Arnold, Kozlov & Neishtadt 
1997, Chapter 5, §2). As in this theorem, the iterative corrections are second-order in the perturbation and are, therefore, 
quadratically convergent. From Hamilton’s equations, the quantity fl = dHpP/dJ is a vector of frequencies and the small 
divisor 1 ■ $7 will occur near a commensurability or resonance. Even when the denominator is bounded away from zero, the 
effect of the denominator is far reaching (see section 2.2.2 below for discussion). 
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This procedure has some favourable features. The iteration produces the series in the new actions and angles directly. Also, 
assuming that the Hamiltonian is available, the computation of the integrands in equation (17) do not require gradients in J; 
that is, the iteration can proceed without computing the frequencies dH/dJ. The updated coefficients are provided recursively 
by (J) t— (J) + (J) and <— Ho + where the initial value of 5 ^^' is simply the first order 

perturbation (see also equation A4). There are, however, several practical problems with the formulation above. The main one 
is that evaluation of equation (17) requires w = w(J,s) which is implicit. This solution be computed straightforwardly using 
the Newton-Raphson method. An alternative approach uses an expansion in the new actions J and old angles w The mixed 
solution avoids an implicit solution of equation (4) to obtain the old angles in terms of the new. A rederivation of the iterative 
solution for the mixed Fourier variables that yields a numerically expedient algorithm for the computation of the expansion 
coefficients using the standard Fast Fourier Transform algorithm without implicit solutions is presented in Appendix A. Each 
choice has advantages depending on the application. 

Finally, the quadratic convergence of the iteration suggests a direct Newton-Raphson evaluation of equation (15). While 
possible, numerical stability seems to require a careful choice of a small but dominating subset of indices 1. Since our goal is 
an automated analysis, we used equation (20) or equation (A 8 ) for the examples in this and in companion papers. 


2.2.2 Discussion 


Mathematically speaking, a term in equation (16) based on the solution from equation (20) with a vanishing denominator will 
result in long-period motion with large amplitude for one degree of freedom. This signals the breakdown of the continuous, 
predictable behaviour of linear perturbation theory. The KAM theory arose as a result of Poincare’s concern about the 
stability of the three-body problem (i.e. the stability of the solar system) resulting from this breakdown. The two systems 
have one huge qualitative difference: dynamical times in the solar system are measured in the billions while dynamical times 
in a galaxy are measured in the hundreds, maybe thousands at most. Owing to our mathematical techniques, the numerical 
characterisation of the perturbation theory failure often requires extremely high-resolution data and time series very long 
compared to the natural frequencies of our dynamical system. This is partly the result of our celestial mechanics bias and 
computational limitations. Galaxies have many orders of magnitude more stars and dark-matter particles than solar system 
bodies of interest. Therefore, even though the available time scales are much shorter, the phase space coverage is so thorough 
that some an observable subset of particles will manifest chaotic dynamics in the time available. This regime begs for study 
and partially motivates this paper. 

Much of KAM theory is concerned with the phase-space domain that admits convergent solutions away from the vanishing 
denominators in equation (20). The theory employs the exponential decay of Fourier coefficients of analytic functions with 
increasing |1| together with the general condition that rational numbers approximate the ratios of Hamiltonian frequencies, 
dH/d\, sufficiently badly that the series converges nearly everywhere for most problems of physical interest under weak 
perturbations. This condition is often written in the following form: 


1 - 


0 tj["l 

L/i-inet 




( 21 ) 


for some (3 > 0 and 7 ^ d where d is the rank of 1. Armed with this formal result, we have confidence that our procedure will 
work much for much of phase space, although we will not explicitly verify the exponential decrease of coefficient amplitudes 
here. However, note that the time-scale mismatch remains; the eternal nature of the torus may be too strong a criterion for 
our galactic systems. We may address this in an intuitive way by noting that the time-scale for a particular term with a 
vanishing denominator will be proportional to the maximum absolute value of the integers in 1 (denoted |il||). Therefore, a 
truncation at modest values of || 1 || combined with the exponential decay of high-order coefficients suggests that the dynamics 
revealed by the algorithm from section 2.2.1 may retain its applicability. We will see examples of this truncation in section 5 
In the standard application of the KAM theorem, one attempts to find a torus subject to the constraint of equation (21). 
This may require some adjustments to the initial conditions in I. Our goal here is somewhat different: we are interested in the 
fraction of tori that break up under the perturbation as a function of phase space volume and therefore, we do not adjust I. 
In other words, we are interested in the converse to the KAM problem: where in phase space do tori not exist? The natural 
implication of failing convergence of the series qi(J) for some J is a broken torus. In practice, failure has two forms: 1) the 
quantity 1 • = e ~ 0 so that the amplitude AWi is too large to continue with the iteration; and 2 ) multiple terms 

of different 1 influence each near one of the separatrices, inducing Smale horseshoes (and chaos thereby) into the dynamics 
around the homoclinic trajectory (as in Holmes & Marsden 1982). Although the KAM theorem does not offer specihc wisdom 
on this, numerical experiments described below and in (alias?) and (alias?) show that failure to converge is most frequently 
an indicator of chaos. 


MNRAS 000, 1-25 (2015) 




8 Weinberg 

2.3 Incremental perturbations 


The KAM procedure is perturbative. Some applications, therefore, may benefit from breaking a large perturbation into small 
increments to better understand the development of irregularity. For example, we might choose to define a perturbation 
strength ej = jelm with j € [l,rn,]. Then, at increment j, the perturbation strength is always of order e/m with 


HoU) 

zz Ho + ^ ^77i, 

m 

( 22 ) 

HW) 

= -77i, 

m 

(23) 

HU) 

= 77o(j)+7fi(i). 

(24) 


Clearly, this could be generalised to a non-linear dependence on e. Then, the algorithm from section (2.2.1) may be applied 
successively to each level j. Explicitly, let the solution at each step j be denoted by VFi(J,w;j), or using equation (2), 
Fi(J,w; j) — J ■ w -f Wi{3,'w\ j). The iteration developed in the previous section results in the following chain of canonical 
transformations: 

[I, w] = [1(0), w(0)] [ 1 ( 1 ), w(l)] [1(2), w(2)] • • • 

[I(m — 2), w(m — 2)] [J-{m — 1), w(m — 1 )] 

[I(m),w(m)] = [J,s]. (25) 

At any step in the chain, the converged perturbed orbit has I(j) = constant with uniformly advancing w(j) at the rate 

n{j) = dH{i{j)-j)/diij). 

For a simple example, consider a chain with a single link, m = 1. We wish to plot the perturbed trajectory for the 
torus with J as a function of time. By construction, the angle variable advance uniformly at the rate fl, so beginning at 
s(0), choose to advance to s{h) = s(0) -I- flh. Using equations (3) and (16), we compute I at s(t = h). Using equations (4) 
and (16), we compute for w at s{t = h). Since the Cartesian values of the phase-space point (x, p) are known since the 
regular orbit is determined by quadrature, we have our phase space point for the torus J at t = fi. For m > 1, the same 
steps are used recursively from beginning at ji = m and moving down the chain by recursion. Equation (4) tells us that the 
w(j — 1) is implicitly defined by [I(j),w(j)]. Solving for w(j — 1), equation (3) gives I{j — 1) explicitly in terms of I(j). 
We then walk down the chain transformations defined by equation (25) for each desired value of [I(j),w(j)] until we reach 
[I, w] = [1(0), w(0)]. The last step provides us with the phase-space point in terms of the unperturbed tori. The penalty of 
using equations (22)-(25) for large m is the large memory needed to store all of the coefficients and intermediate tables of 
actions and angles at each j. See (alias?) for additional discussion of this topic and its application. 


2.4 Numerical considerations 
2.4-1 Numerical partial derivatives 

The computation of new perturbed action-angle values along with the computation of the new frequencies requires action-space 
partial derivatives of the expansion coefficients from equations (17) and (20) using equations (3), (4) and (5). These numerical 
computations rely on either finite-difference or functional approximation techniques, and as usual, the best choice is problem 
dependent. We will use the notation defined in Tremaine & Weinberg (1984) for action variables in a spherically-symmetric 
gravitational potential: 7i and I 2 denote the radial and azimuthal action in the orbital plane, and /s denotes the projection 
of I\ on the axis perpendicular to the equatorial plane (i.e. the 2 axis). Depending on the perturbation, the perturbed 
actions Ji, J 2 , J 3 will no longer have these clear interpretations. Eor a perturbation to a thin disk by a vertically symmetric 
disturbance, we can restrict our consideration to the action pair (/i,/ 2 ) and consider only a two-dimensional distribution of 
action points. In general, an initially rectangular grid of nodes (/i,/ 2 ) will not rectangular after the perturbation is applied. 
Therefore, we can not rely on the standard algorithms that assume a rectilinear grid for (Ji, J 2 ) for interpolation. Eor this 
paper, we explored two techniques to handle the expected irregular distribution of nodes (Ji, J 2 ). Eirst, the modified Shepard 
algorithm (Shepard 1968) uses a Voronoi tessellation to provide a metric on the unstructured grid and inversely weights the 
contribution of nearby points by distance within a specified radius. Alternatively, one may use a k-d tree to identify points 
within the specified radius using a nearest-neighbour search algorithm. The latter is slightly more efficient in greater than 
two dimensions. Secondly, Riachy, Mboup & Richard (2011) find an optimal set of coefficients in the least-squares sense to a 
multidimensional Jacobi polynomial approximation in an irregular region. The Shepard algorithm performed more robustly 
in tests and we use this approach for this and the companion papers. 

The third action, /a, is conserved for an axisymmetric but otherwise three-dimensional perturbation (such the effect 
of a flattened perturbation on spherical bulge and halo). Therefore, we may consider slices of constant /a individually and 
break the problem into two pieces. We first define planes of constant 7a = 7/z. The maximum value of Is may be defined by 
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choosing an effective circular radius: Vdrc and setting I 3 = rdrcVcirc- Since I 2 ^ I 3 and I 3 = 72 cos(/3) where (5 describes the 
inclination of the orbital plane relative to the z = 0 plane, once the value of I 3 is fixed, the values of I 2 > I 3 determine /3. In 
other words, the orbits defined on the I 1 -I 2 plane at fixed I 3 have larger inclination off the disk plane as I 2 increases. This 
is a slightly odd way of representing phase space but is convenient for numerical derivatives. For example, fla = dH{I)/dl 3 
may be computed by interpolating H for the desired values of (7i, I 2 ) and using finite-difference derivative between I 3 planes 
to compute 977(1) 79 / 3 . Because I 2 ^ |73|, the computation of derivatives at the lower boundary in I 2 may require special 
treatment such as extrapolation. 


2.4-2 Computing the action and angle transforms 

We truncate the action-angle series by choosing a maximum index in each dimension, These indices should be chosen 

sufhciently large to obtain an accurate representation of the orbit by the trigonometric polynomial. For eccentric orbits, 
whose expansions include terms with large ||1|| owing to the high instantaneous frequencies at pericentre, the required value 
of lj,max may be much larger than the desired number of retained commensurabilities. However, overly large values of lj,max 
are computationally costly, both in CPU time and memory usage. For computations here, we find that Ij^max ^ 16 is a good 
compromise. We use the efficient FFTW package^ for the computing the discrete Fourier transform using the fast-Fourier 
transform algorithm. 

One of the pair of transform equations (equations 3 and 4) is implicit, depending on whether we use the angle expansion 
in new variables (as in equation 16) or in old angle variables (as described in Appendix A, equation Al). As described at 
the close of section 2.2.1, the mixed variable form that follows the Type 2 generating function explicitly is computationally 
more efficient. In this case, transforming from the new action to the old action given the old angle is a simple application of 
equation (3) given uji from equation (A7) or equation (A8): 

1 = J+^ = J-bi^lt:u,(J)e‘‘"'. (26) 


The angle transform is more complicated in two ways. First, equation (4) requires an action derivative of uJi(J), which we 
perform using the modified Shepard algorithm as described in section 2.4.1. Second, equation (4) is now an implicit expression 
for s in terms of w: s = s(J,w). Fortunately, this is efficiently solved using the Newton-Raphson method as follows. Define: 



Q(w) = w-1-^-^e 

1 

(27) 

and 

1 

(28) 

Then, the Newton-Raphson iteration for 

w becomes 



Wj = Wj_i -1- M~^ ■ ( ® “ Q(wj-i) J . 

(29) 


This converges to machine precision in under ten iterations in practice. 


2 . 4.3 Time dependence 


The formalism of section 2.2 may easily be extended to time-dependent perturbations. Similar to the transformation of the HJ 
equation to an algebraic equation using action-angle transforms, we may treat the time variable using a Laplace transform. 
The special case of a simple harmonic time dependence, such as that resulting from a constant pattern speed with frequency 
mflp, a Fourier transform is sufficient. The resulting transform has a single term and the denominators in equation (20) 
become: 


9rrN 977^"' 

wiinew I 

“9J ^ 


— mflp. 


(30) 


We sum over all driving frequencies j if there are a finite number of discrete constant drivers mjQpj. If the time-dependence 
is aperiodic, the resulting transform results in a continuous frequency-like Laplace transform variable s: 


1 


9777 


9J 


1 


9777 


9J 


+ is. 


(31) 


The inverse Laplace transform required for the fully general time-dependent Hamiltonian requires an additional grid of 
quadratures but otherwise, the algorithm follows as in section 2.2. In the example below and in (alias?) and (alias?), we 
will assume at most a single driving frequency. 


^ http://www.fftw.org/ 
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3 SOFTWARE IMPLEMENTATION 
3.1 Tori computation 

The ideas from the previous sections have been implemented in an object-oriented C-|--|- code. Readers with little interest 
or experience in object-oriented programming (OOP) may skip this section without consequence. In short, an object or class 
is a data structure and a set of functions. OOP imbues the relationship between be classes with the same relationships and 
interrelationships that their counterparts have in the physical or mathematical world. Thus, OOP provides a framework for the 
computational scientist to build and extend their numerical machinery (see Pitt-Francis & Whiteley 2012 for an illustrative 
introduction). Because OOP enforces the natural relationships between objects, it also helps prevent mistakes. The specific 
OOP design features for the code that implements the algorithm in section 2.2.1 are as follows: 

(i) The Hamiltonian torus concept is embodied by a single class whose instances track the action vector, the coefficients 
of the Wi expansion, and the derivatives of Wi necessary to compute the angles at any desired time and at any level of the 
perturbation hierarchy. 

(ii) Time-dependent perturbations will require multiple harmonic orders or frequency subspaces. These may be used to 
track periodic or aperiodic temporal perturbations using Fourier or Laplace, respectively. These subspace domains are fixed 
by the definition of the perturbation. The torus class instances keep lists of values indexed by each subspace. 

(iii) A container class keeps track of the full collection of action-vector nodes and associated tori for each perturbation. 
This container implements the interpolation and differentiation of the possibly irregularly sampled action points that enables 
the computation of the perturbed trajectories for each new torus. 

(iv) This container class may have any number of derived classes that may initialise the individual nodes with action values 
based on specific symmetries or special geometries in phase space. Since the calculation begins with regular orbits, a particular 
derived class implements the computation of the analytic mapping between Cartesian and angle-angle phase-space vectors for 
each geometry of interest. The astronomical applications here and in the companion papers assume axisymmetric or spherical 
unperturbed models. 

(v) The container class of nodes has been designed to handle a three-dimensional unstructured grids, two-dimensional grids, 
and three-dimensional grids restricted to two-dimensional slices with fixed I 3 . During development, we tested the general three- 
dimensional grid using a computation that was restricted to two-dimensional slices. As expected, the implementation of a 
full three-dimensional grid is the most expensive and noisiest, requiring an order of magnitude more tori to achieve the same 
accuracy than the equivalent restriction to two-dimensional slices. 

(vi) The container class keeps track of and attempts to diagnose any problems in the interpolations as the computation 
proceeds. Specifically, for strong perturbations, the induced oscillations for one zui term on another may drive excursions 
beyond the initial action grid. Most often, these happen near phase-space boundaries for diverging iterations typical of broken 
tori. 

(vii) Finally, these computations often represent an investment in both human and computing time. All intermediate results 
are stored to full precision for later reuse using the BOOST (http://www.boost.org) serialisation library. 

For each newly-defined computation, we begin by choosing an appropriate value of Ij^max for each dimension in the 
subspace. This finite set of “quantum numbers” results in a cubical lattice of For two-dimensional problems, the 

lattice is a square. We have explored the dependence of the coefficients on the selected value by increasing Imax by factors of 
two until convergence reached. For problems explored to date Imax ^ 16 has been sufficient. Also, significantly larger values 
require special accommodations memory requirements. 


3.2 Computation of Lyapunov exponents 

The Lyapunov exponents, Xj for j = 0, ..., n — 1 where n is phase-space dimensionality, measure the exponential divergence 
rate of initially nearby phase-space points. The so-called standard algorithm (Wolf et al. 1985a) exploits the multiplicative 
ergodic theorem (Oseledets 1968). The import of this theorem may be motivated by the following example. Consider the 
mapping defined by the continuous solution of the linearized equations of motion, u = G(x, t) • u where Gij = dFi/dxj. This 
gives us a tangent-space mapping from t — 0 to t = h: u{h) = T(i[x(0)]u(0). If our tangent map Th were constant for all time, 
it would be a trivial matter to compute the eigenfunctions and eigenvalues once, and compute the subsequent evolution at any 
desired time. We cannot do this because T is time dependent, and therefore we would have a new eigenvalue problem at every 
time t. Oseledets’ theorem tells us that the generic case behaves as if T were constant and that the Lyapunov exponents are 
related to the eigenvalues of the the mapping matrix: Xj = limt_^oo t~^ log(|T • Uj|/|ujj) for appropriately chosen eigenvectors 
W- 

The major shortcomings of the algorithms based on this theorem are that the accuracy of the estimates depend sensitively 
on the numerical method. Because of this, there is a very large literature describing improved Lyapunov exponent algorithms. 
However, for our purposes here, accurate exponents are not required. We will use only the standard algorithm to discriminate 
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between exponents that exceed some positive threshold value and those that do not. Also, following common practice, we will 
use the continuous Gram-Schmidt normalisation procedure (Christiansen & Rugh 1997) to estimate the various exponential 
subspaces to prevent the computational singularity that results from the direct application of the Oseledets theorem. There 
are many newer approaches with virtues for different applications (see Darriba et al. 2012, for a review); the choice of the 
standard algorithm for this application may not be optimal. 

Our ODE solver is the standard RK4 stepper. We experimented with a 6th order symplectic Runge-Kutta (using the 
Yoshida 1990 trick) and the Runge Kutta Felhberg algorithm (e.g. Cash 1985), and the adaptive Bulirsh-Stoer (B-S) method 
with a specified relative error (e.g. Press et al. 1992). The non-sympletic B-S scheme is much more efficient owing to the 
multiple length scales of the halo/spheroid and disk. Relative energy changes of one part in 10® and relative changes of 
one part in 10® are obtained with time maximum time steps of 0.0001. However, an externally chosen fixed time-step using 
the standard RK4 yield the most rapid convergence of the Lyapunov exponents on average. 

The RK4 time-step, 5T, is chosen as follows. We begin with a time step selected using the typical n-body simulation 
algorithms: a fraction of the characteristic dynamical time determined by either the logarithmic derivative of the position 
or velocity. Then, we solve the equations of motion numerically to find the smallest ratio of velocity to force. We use 1/100 
of this value to integrate the equations of motion and its tangent space to determine Lyapunov exponents. This fraction is 
both increased and decreased to check for convergence of the derived exponents. Typically, this procedure results in a time 
step one order of magnitude smaller than that chosen by the traditional algorithm used for n-body simulations. In addition 
to the ODE-solver time-step, one must choose the number of individual time-steps between Gram-Schmidt orthogonalisation 
of the principle directions. After some experimentation, we choose this number to be 1/3 of the mean orbital time. Our 
implementation was tested successfully using the Lorenz equations which has well-known Lyapunov values (Wolf et al. 1985b). 

Convergence of the maximal Lyapunov exponent. A, for small value of |A| can be logarithmically slow. That is, if the true 
value is A = 0, the most positive exponent approaches zero proportional to 1/N where N is the number of steps 5T. For the 
work here, we are less interested in the value of A but whether A > 1/T for some characteristic time T. 

Our algorithm is as follows: 

(i) Discard the first NiniUai of the recorded Lyapunov values. Typically, we choose a series j G Ntotai = 2 x 10® steps of 
size 5T discard the first half of the series; that is, we set NiniUai = Ntotai/2. 

(ii) Order the exponents Xoj > Xij > • ■ • > Asj. For Ao, compute the linear regression of log Aoi with log/ to estimate 
the slope m. 

(iii) Perfect convergence to zero yields \m\ = constant = 0(1) while perfect to a non-zero positive value yields m = 0. If 
the computed slope m is greater that —1/2 we assume that the most positive exponent is zero. 

The phase space vector is computed from the background potential for each torus in the grid (as described in section 
3.1). In this paper, we are typically analysing of order 10'* different initial conditions for multiple numerical experiments, so 
examination of each by eye is not feasible. However, we have done this for some test cases and failure modes in this algorithm 
warrant some discussion because they appear to indicate interesting dynamical situations. First, we noticed that the run of 
individual Afc with step j do not maintain their initial ordering. For example, it is not unusual for the value of Afe to converge 
a different rates. Secondly, and more importantly, individual orbits that appear to be converging to positive value Ao = e > 0 
at small j will begin to converge to zero logarithmically at large j. Occasionally, the opposite behaviour is seen: a value of 
Ao oc 1// at small j may suddenly increase and converge to decreasing value of Ao = e > 0 at large j. 

After checking for and not finding any implementation errors, we concluded that this behaviour is consistent with the 
dynamics under study. In particular, an orbit in a chaotic region will exhibit exponential sensitivity. However, the apparent 
diffusion in action space may result in an orbit “sticking” in a region of apparent regularity for some time, resulting in a 
decreasing positive value of Ao. This well-known behaviour is the result of an initially diffusing orbit encountering regular 
islands around separatricies and at the boundary of stochastic sheaths. Leaving a detailed study for later, we note that this 
tendency towards stickiness in stochastic layers surrounding separatricies of strong resonances undermines the utility of using 
Lyapunov exponents to characterise the dynamics in these stochastic layers. Also note that these same orbits that exhibit 
non-uniform convergence often coincide with values of broken tori identified using the numerical KAM method, suggesting 
that these are regions of weak chaos (i.e. sub-exponential sensitivity) that occur near strong resonances. Phenomenologically, 
these weakly-chaotic trajectories appear regular throughout much of their orbital phase but may switch morphology suddenly 
near heteroclinic points. So, although they may appear to have Lyapunov exponents close to zero, they give rise to orbits 
with varying behaviour. 


4 INTERPRETATION OF BROKEN TORI 

Hamilton’s equations for a Fourier-series perturbation with a single dominant term (i.e. a single ■uj\ in equation 16 or Al) 
may be solved directly without difficulty. However, difficulties do arise when several terms contribute comparable power to 
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equation (16). For example, consider two such terms {A and B) with 1 a 7 ^ Is- As described in section 1.2, each term alone is 
analogous to a single pendulum with a single degree of freedom. Consider an initial condition near the unstable equilibrium 
of one of the pendula. The initial conditions must be finely tuned for the trajectory to ‘climb’ to the unstable equilibrium. As 
the trajectory approaches this equilibrium, the force of the primary pendulum vanishes, and then, the other pendulum can 
influence the trajectory, causing apparent stochastic behaviour. Depending on the state of the full system, this influence can 
cause the trajectory to enter the either the positive or negative rotation zones or the libration zone of the primary pendulum. 
In summary, the two interacting pendula have resulted in the destruction the torus described by the action conjugate to the 
angular degree of freedom describing the unstable equilibrium. In some cases, for very strong perturbations, resonances from 
harmonics of the primary resonance can interact. This dynamics leads to the resonance overlap criterion (Chirikov 1979). 
This will be familiar in solar-system context as orbit instability induced by the overlap of mean-motion resonances. 

The consequences of this broken torus depends on the problem at hand. If one of the pendula corresponds to a primary 
resonance, noise near the resonance can cause the ensemble of trajectories representing the entire torus (e.g. initial conditions 
with the same initial actions but uniformly spaced phase values) to change their actions depending on their phase. Moreover, 
the problem can be complicated by many Fourier terms simultaneously interacting and generating wide regions in phase 
space that exhibit exponentially sensitive stochastic behaviour. The erratic behaviour can be dramatically exacerbated if 
the unstable trajectories generated by several Fourier terms coincide in phase space. In this case, trajectories may move 
from resonance to resonance. In some cases, the influence of individual components may be too plentiful to discriminate. As 
previously described (see section 3.2), the degree of stochasticity is most often characterised by the product of Lyapunov 
exponent, A and the characteristic dynamical time, r. Values of C = Ar ~ 1 are called strong chaos and values of C < 1 are 
called weak chaos (e.g. Zaslavskii et al. 1991). 

The method highlighted in section 2 is based on a perturbative solution including all the lower-order terms 1 within some 
lattice with cardinality lj,max- We assume that failure to find a solution leading to a new torus implies a broken torus: the 
absence of an analytic predictive description of future time evolution of an ensemble of nearby phase-space points. Moreover, 
although broken tori and Lyapunov exponents are related, an exploration of phase space using each method reveals different 
information about the underlying dynamical structure. 

For example, loci of broken tori in a phase space perturbed away from regularity may result from a specific strong 
resonance. Orbits of different morphology reside on either side of or in islands surrounding this resonance. With increasing 
perturbation strength, these loci broaden as additional terms conspire to break tori. However, although a broad region of 
broken tori implies regions of positive Lyapunov exponents, this does not imply that C > 1, and conversely, large C may occur in 
narrow regions around separatricies, as in the standard map. As we will see below and in (alias?) and (alias?), regions of high 
C often appear in regions of broken tori but not exclusively. In some cases, the largest values of C occur between the widening 
loci of primary resonances. The loci of original primary resonances sometimes remains a local minimum in C, presumably 
owing to regular islands spawned by these resonances, as suggested by the morphology of Poincare surface-of-section plots for 
the standard Hamiltonian mapping problem. 

4.1 A simple dynamical model 

In this section, we investigate the evolution of a two-pendulum model with two non-zero terms in equation (16) to explicitly 
contrast the one-degree of freedom problem that underpins secular evolution with the multiple degree of freedom problem 
that admits chaos. For an explicit example, let us consider the Hamiltonian 

7L(I,w) = H„(I)-bti 7 A(I)e”‘-^-"'-ttus(I)e"«'"', (32) 

where (I,w) are the action-angle variables for the unperturbed Hamiltonian Ho{T) with I = (/i,/ 2 ) and w = (wi,W 2 ). The 
restriction to a Fourier single term in the expansion series results in classical secular perturbation theory. Occasionally, the 
linear solution described above may diverge owing to the resonant denominators, but a non-linear solution can be found 
most often for a single term. However, when one or more additional Fourier terms are added to the Hamiltonian, irregularity 
arises. The importance of the irregularity to the solution in its astronomical context depends on the relative magnitude of the 
expansion coefficients and the order of commensurability. We will quantify the order of resonance as |{lj{ = suplZ^}. 

It is well-known that chaos arises in a stochastic layer near the homoclinic trajectory (e.g. Reichl 2004). So, we focus on 
the solution of the perturbed Hamiltonian near the homoclinic trajectory for a single perturbation term, say wa 7 ^ 0, wb = 0. 
Motivated by orbit averaging, we transform to slow and fast variables defined by the following Type 2 canonical transformation: 

F={U-^)Is+wJf. (33) 

Near the commensurability defined by 1a ■ w = 0, equation (4) implies that the new angle Ws = 1a • w will change slowly 
and the new Wf = Wi will change rapidly with time. The expressions for B and If follow similarly from equation (3). The 
choice of the second term has a fair bit of flexibility in this application; the most important constraint is that \wf\ >0 when 
I tils I « 0. These conditions can be made rigorous (Lochak & Meunier 1988) but are not needed for this example. 
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Figure 1. Poincare surface-of-section plot with zua = 1 and ujb = 0 showing action Is versus angle Ws for Wf = 0 for approximately 
10'^ recurrences. The action of the homoclinic trajectory is Is = ±2. Trajectories with |7s| < 2 and \Is\ >2 are in libration and rotation, 
respectively. 


To define an explicit form for the Hamiltonian in equation (32), we expand -Ho (I) to second-order around a set of 
unperturbed initial actions (/so,-^/o)- Ignoring the cross term and completing the squares, we obtain: 

H{J, w) = C(Io) + <^^+ (34) 

where C(-) is an uninteresting constant and a and (3 represent d^H/dJ^ and d^H/dJ^ respectively, the definitions of Js, Jf 
include the offsets from initial conditions and completing the squares, and 7 and 5 are Ib 2 /Ia 2 and Ibi—Ib21ai/Ia2 , respectively. 

For a simple example, we begin with initial conditions near the homoclinic trajectory defined by a = 1, /? = 1/2, 
VJA = 1,^B ~ 0 with 1 a = (—1,1) and Is = (1,1). This implies that 7 = 1 and 5 = 2. The tub = 0 surface-of-section 
(SOS) plot (Fig. 1) reveals that standard phase space of the regular, integrable pendulum. For zub > 0, the homoclinic 
trajectory broadens owing the action of the non-resonant term on on the resonant term near the unstable point. The SOS 
plot for zjb = 0.01 and tub = 0.1 are shown in Figs. 2 and 3, respectively. Typical values of vobIz^a for astronomical 
perturbations described below and in the companion papers are between 0.1 and 0.5 so Fig. 3 is the more typical situation. 
For values ■ujbI'zoa = 0(0.1), much of the libration is chaotic with a relatively wide chaotic zone around the original homoclinic 
trajectory. The dynamics of models like this are nicely described in chapters 2 and 3 of Zaslavsky (2007). 

The model of equation (34) is today example only; most realistic expansion series have approximately 5 terms with 
amplitude ratios with the primary between 0.1 and 1.0. Many may interact simultaneously leading to complex chaotic dynamics 
near the primary resonant orbit. This simple model helps motivate the prevalence of chaotic behaviour in the astronomical 
examples in the companion papers. Additional investigation of the dynamics of this simple model reveals a wide variety of 
features found in the standard Hamiltonian mapping, 

-20„=-ifsine^ (35) 

with A" > 1. In particular, one finds a stochastic ‘sea’ with a wide variety of islands whose location change quickly with 
small changes in action within the stochastic layer. This leads to the behaviour of Hamiltonian intermittency: the appearance 
of both stochastic and regular behaviour in the same trajectory. This results from trajectories being held up at the border 
of stable islands, sometimes called “stickiness” (i.e. Lee 1988; Lai 1999). These dynamics may lead to heavy-tailed escape 
distributions similar to those caused by the first-order Fermi effect (Zaslavskii et al. 1991). The relevance of possibly pervasive 
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Figure 2. As in Fig. 1 with rog = 0.01. Although the amplitude of the second term is only 1% of the primary term, a region around 
the homoclinic trajectory is chaotic. Also note the distorted resonant orbit for Ig = 0.4. 


weak chaos in shaping the galactic velocity distributions and its relevance to so-called radial migration motivates additional 
investigation. 


5 EXAMPLE: CIRCULAR RESTRICTED THREE-BODY PROBLEM 
5.1 Introduction 

The three-body problem is one of the oldest and best studied in dynamics, although surprises continue to emerge. Here, we 
will consider application of the numerical KAM method to the restricted three-body problem. The main goal for this example 
is the comparison our diagnosis of broken tori using the nKAM method with the characterisation of chaos using Lyapunov 
exponents for a familiar problem. 

We take units GMq = 1, consider the perturbation to be a circularly orbiting Jupiter-like planet with Mp = 0.001 and 
a = 1, and explore the motion of a coplanar test particle. We examined two types of perturbations. In the first, we assumed 
that the unperturbed potential is that of a central point mass fixed at the origin (i.e. the Sun) and the perturbation is an 
orbiting point mass (i.e. Jupiter). In the second, we assumed the perturbation is the difference between the unperturbed 
potential and that generated by the full two-body solution including barycentric motion. The results were nearly the same 
for both perturbations for harmonic orders m = 1 through m = 3, so we will use the fixed barycentre perturbation for the 
remainder of the calculations. 


5.2 Tori computation 

Here, we describe the results of the torus construction outlined in section 2. Fig. 4 shows a plot of broken (coloured circles 
and squares) and unbroken (blue dots) tori as a function of radial and azimuthal actions, li and I 2 , for the Jupiter-like model. 
The original grid of radial and azimuthal actions is chosen by choosing a range of minimum and maximum energy defined by 
the circular orbits between rmin = 0.05 and Vmax = 120.0. Let E{rc) and Jmax{E{rc)) be the energy and angular momentum 
of the circular orbit at r^. For each energy value in the range [E{rmin), E{rm.ax)], we compute the minimum and maximum 
values of Ji and I 2 such that 0.001 < I 2 /Jmax{E) < 0.999. Phase-space points outside the envelope of this region are shown as 
white. The legend has been restricted to low-order resonances for clarity. Other broken tori, for higher-order tori, are shown 
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Figure 3. As in Fig. 1 with rog = 0.1. The amplitude of the second term is only 10% of the primary term; most of the libration zone is 
filled with chaotic trajectories, although regular initial conditions can be found still. For example, trajectories with = 2.0, —1.9 near 
the original homoclinic trajectory were chosen to display embedded islands. An example of a resonant island trajectory is Is = —2.5. 


as black points. To help interpret the results, note that the two-body (Kepler) problem in the adopted units has the following 
simple representation as a function of I: the energy is 


H = 


E = - 


1 1 
2 (7i -b hy ’ 


(36) 


the eccentricity is 


€ = 



(37) 


and the semi-major axis is: 


a = { Ii + I2 


(38) 


Therefore, diagonal lines at constant 7i -|- I 2 are lines of constant energy and semi-major axis. Loci of constant eccentricity 
are the lines Ii = [(1 — — l]l 2 \ the diagonal I\ = I 2 has eccentricity e « 0.866. Because the frequencies are independent 

of e at fixed E, these diagonals are also lines of constant orbital frequency and therefore loci of resonances as well. Initial 
phase-space values were selected from a grid of 7i and I 2 values chosen between upper and lower energy bounds. The finite 
steps in each action value yields non-rectangular boundaries in the I 1 -I 2 plane. 

The method described in section 2.2 is most efficient when the Fourier expansion series converges quickly, that is, when 
the value of Imax is small. However, for this planetary perturbation, orbits with increasingly large e can pass close to the 
perturbing planet for increasingly large values of I\ (or energy E). Convergence of the Wi series would require increasing large 
values of Imax- Clearly for modest fixed values of Imax, this method will only produce accurate tori estimates for I 2 ^ Ii- 
Nonetheless, our example below will show that we can still obtain interesting insight into the dynamics with a small finite 
cube of 1 pairs. 

A number of features are immediately evident. First, the broken tori cluster around loci corresponding to principle low- 
order resonances. We identify the coefficient Wi with the largest magnitude as the principle resonance of each torus; the 
low-order principle resonance “quantum” numbers (lijh) are indicated in the legend. The dominating terms in the series 
should reflect the underlying drivers of the perturbation, however, the identification of principle resonances is solely based 
on amplitude in Figs. 4 and 5. In many cases, the lowest-order term is only slightly smaller than the largest term with 
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Figure 4. The broken and unbroken tori using the nKAM method for |Zi| ^ 8, IZ2I ^8 and azimuthal harmonic orders Z3 = m = 0,..., 3 
in the action plane. The legend indicates primary components Wi in the form {li,l2) for |Zj| < 3, broken tori of higher order (|/j| ^ 3) 
are indicated by black dots and unbroken tori are indicated by blue dots. 


higher-order 1. Therefore, the identifications in this figure may not precisely the morphological character of the broken torus 
but still provide some clues to the important features of the dynamics. Two or more terms Wi of comparable amplitude are 
most likely necessary for stochasticity. If desired, one may further study the details of specific broken tori using the direct 
integration approach described in section 4.1. 

The tori with a < 1 unbroken except at high eccentricity, e > 0.95. For reference, the locus I 2 > 47i have e < 0.6 and 
I 2 = Ii has £ = 0.866 (from equation 37). The swath of broken tori with Ji > 0.5, I 2 < 0.3 have been spot checked by direct 
integration and many of these do indeed exhibit unusual stochastic behaviour “by eye” in the form of eccentricity changes 
and sometimes irregular precession. The primary resonance identification suggests that resonance spreads out in to the cloud, 
presumably aided by other high-amplitude subdominant terms Wi. The locus of broken tori in the upper right of the diagram 
are induced by the 3:1 and 4:3 resonances and other mean-motion harmonics with other high-order secondary resonances. 
The SOS in planet’s rotating from for the upper-left most torus in this locus is shown in Fig. 6. 

For comparison. Fig. 5 shows a grid of tori at lower action-space resolution but larger maximum energy and radius. This 
figure shows that the damage to the tori by the Jupiter-like planet is localised to the vicinity of the planet and eccentric planet 
crossing orbits, perhaps as expected. The upper locus is, again, the result of the 3:1 and 4:3 primary resonances. Broken tori 
in this locus are affected by significant stochastic behaviour as evidenced by direct integration. 

Direct integration broken tori at small Ii and 0.5 < /2 < 1-7 show less dramatic stochastic behaviour. There are a large 
number of distinct groups of broken tori in Fig. 4 that might be of interest for various dynamical situations. Most of these 
are likely to have been investigated previously by planetary dynamicists. Our main goal, here, is to make contact with a 
well-established problem. Key to this, is understanding whether the failure of the numerical KAM procedure to produce a 
torus implies an irregular, stochastic orbit. To better assess stochasticity, we turn to Lyapunov exponent computation using 
the methods described in section 3.2. 

For direct confirmation of broken tori, the Melnikov method has become the standard tool for detecting splitting of 
invariant manifolds for systems of ordinary differential equations close to integrable ones. The Melnikov method for two 
degree of freedom Hamiltonians described by Holmes & Marsden (1982). Their method estimates the magnitude of the chaos 
by measuring the fraction of the energy surface where chaotic motion occurs. This idea will be explored in a future contribution. 
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Figure 5. As in Fig. 4 but for an action grid with larger maximum energy (larger semi-major axes), showing that broken tori are 
localised to the vicinity of the planet. 



Figure 6. SOS plot for broken torus in upper-left most corner of the 3:1 resonance locus in Fig. 4 with {7i,/2) = (1.132,0.837). 
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Figure 7. Density of action points from Fig. 4 with positive Lyaponov exponents determined using the algorithm described in section 
3.2. The data is somewhat under-smoothed to better reveal the fine-scale structure. 



Figure 8. As in Fig. 7 using the action points from Fig. 5. 
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Figure 9. A chaotic orbit with 7i = 0.185 and I 2 = 0.742 (red). See Fig. 4 to locate this broken torus. The orbit of the Jupiter-like planet 
is also shown (green). The primary broken torus has commensurability Zi = 1, Z 2 = 0, la = 1, with secondary terms Zi = 1, Z 2 = —1, Z 3 = 2, 
and Zi = 1 , Z 2 = 1 ,13 = 0 . 


5.3 Lyapunov comparison 

We use initial conditions generated from the unperturbed tori in Fig. 4 to estimate positive Lyapunov exponents using the 
method outlined in section 3.2. To recall, the Lyapunov exponents for each orbit are computed using the tangent-space method 
described by Christiansen & Rugh (1997). To diagnose convergence, the time series of exponent estimates, computed from a 
running finite-time average, are analysed using linear regression computation. If the largest exponent approaches zero with a 
slope m > —1/2, the exponent is assumed to be zero; otherwise, it assumed to be positive. The relative density of positive 
Lyapunov exponents indicating stochasticity are shown in Fig. 7. A value of 1 (0) indicates that all (no) neighbouring orbits 
have positive Lyapunov exponent. This figure is computed by smoothing the irregular distribution of initial (/i,/ 2 ) values 
onto a regular grid with spacing 0.01 using a Gaussian kernel with width 0.01 in action value. The initial smoothing value was 
estimated by a cross-validation analysis (Silverman 1986) and manually adjusted downward by a factor of 5 to better reveal 
some of the intrinsic structure common to Fig. 4. 

The comparison of Figs. 4 and 7 reveals both similarities and differences. Overall, the regions of broken tori coincide 
with regions of positive Lyapunov exponents. The 3:1 resonance locus at 7i « /a « 1.8 also appears as a region of strong 
stochasticity. Not present in Fig. 4 are the horizontal loci of high-density positive Lyapunov exponent with 7i « 1.4 and 
I 2 > 0.5. Additional investigation revealed these are caused by high-order secondary resonances not included in the Wi series. 
These are eccentric trajectories that come close to the planet at pericentre; their inherent stochasticity would only be revealed 
by the nKAM method for large ||1|{. Similar a horizontal locus with I 2 < 0.5 is present in Fig. 4. The direct integration of 
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the equations of motion for a particular initial condition in the broad region of broken tori is shown in Fig. 9; changes in the 
argument of periapsis and eccentricity results from close encounters with the perturbing planet, as expected. 

Figs. 5 and 8 compare the same model with lower action resolution but larger maximum energy, large radii orbits. This 
further illustrates the major features of the findings here. First, the low-order resonances give rise to stochasticity spreading 
from these resonances in the diagonal region described by 0.5 ^ /i + ^2 ^ 1-5. Second, these stochastic zones spread for high 
eccentricity orbits (7i > I 2 ) forming a horizontal band of chaotic orbits with I 2 ~ 0.3. This region has eccentricity e > 0.9. 
This band results from low-order terms interacting with high-order terms. Fig. 8 reveals an additional band at larger I 2 , 
‘parented’ by higher orders resonances. 

The locus of broken tori parented by the 3:1, 4:3 and other harmonics primary appears to demarcate the the upper 
horizontal tongue of high-eccentricity (7i 3> I 2 ) positive Lyapunov exponents. As in the case of Fig. 7, the lower horizontal 
loci have no detected counterpart in the numerical KAM results, presumably owing to the truncation in |jl|| as previously 
discussed. The region of positive Lyapunov exponents is bounded from above by the locus of eccentric orbits with pericentre 
values of ~ 3 and bound from below by the locus with apocentre values of ~ 0.5. 

Overall, we find that the nKAM method and the Lyapunov exponent analysis provide a consistent picture of the broken 
tori in the restricted three-body problem for a Jupiter-like planet. The nKAM method provides useful dynamical details, 
such as which resonances are responsible for the stochastic behaviour around the separatrix. This is done by recording the 
dominant terms in the series defining the perturbation (equation 16). These terms may be later used later for detailed analyses, 
as desired. The Lyapunov exponent analysis gives a broad-brushed view of stochasticity and is not restricted to low-order 
resonances. 


6 DISCUSSION AND SUMMARY 
6.1 Main results 

This paper describes a method for estimating the existence and destruction of regular orbits or tori under a perturbation. All 
orbits are assumed to be regular before applying a perturbation. This is most easily realised by choosing a standard model in 
Stackel form (Stackel 1892; Eisenhart 1934; Binney & Tremaine 2008); spherical and two-dimensional axisymmetric potentials 
are trivial examples. Classical Hamiltonian perturbation theory isolates the response of individual resonances by identifying 
the particular slowing varying degree of freedom identified with the commensurable frequencies at resonance and averages over 
the rest. However, this procedure fails to include important dynamics near the resonance (i.e. close to the separatrix of the 
principle resonance) owing to non-negligible perturbations by the averaged terms that have been cast away by the averaging. 
The new method proposed here uses the Hamilton-Jacobi equation to extend the reach of canonical perturbation theory by 
including the mutual interactions of excited resonances. It is now well known that these terms lead to dynamical instability 
through classical indeterminism or chaos (e.g. Zaslavsky 2007). By retaining many terms in the original series rather than 
dispatching them using the averaging theorem, we develop a perturbation theory that includes chaos. 


6.2 Motivation and interpretation 

The new method is motivated by mathematical machinery used in some proofs of KAM theory (e.g. Poschel 2001, section 
3) and their correspondence to techniques used by this author previously. Specifically, both the KAM procedure and the 
standard secular perturbation theory attempts to solve the perturbed Hamiltonian equation using a Fourier analysis. The 
KAM procedure defines an iterative procedure to obtain a new set of action-angle variables as the solution to a perturbed 
Hamilton-Jacobi equation. The standard secular perturbation theory uses the averaging theorem to derive one-dimensional 
equations of motion for each resonance of interest. The essence of the method described in section 2 is a quadratically 
converging iteration for the coefficients of a series approximation to a canonical generating function that defines a new set 
of action-angle variables under the perturbation for all terms simultaneously. That is, the converged coefficients give us a 
canonical generating function that solves the perturbed Hamilton-Jacobi equation and provides new action-angle variables 
(i.e. new tori). Moreover, the KAM theorem suggests that it is sufficient to retain some finite subset of coefficients for some 
modest perturbation strength. This justifies in part the truncation of the expansion series that is necessary for practical 
numerical computation. 

When the coefficients for a particular value of the action vector, I, fail to converge, the cause is rapid oscillation near 
the unstable equilibrium of the primary resonance induced by one or more secondary resonances. This is an example of the 
so-called homoclinic tangle discovered by Poincare. Specifically, the perturbation splits the stable and unstable parts of the 
infinite-period resonant trajectory; the frequency of their mutual intersections intersect increase rapidly near the original 
unstable points of equilibrium. Wiggins (1992) describes the dynamics of these tangles and their consequences in exquisite 
detail. From a heuristic perspective, the KAM theorem tells us that the onset of chaos is a result of coarse graining the 
dynamics generated by the homoclinic tangles. The consequences of these tangles is shown in the stochastic layer that appears 
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around the unperturbed homoclinic trajectory (Fig. 1) as the amplitude of the perturbation is increased (Figs. 2 and 3) in 
our toy model from section 4.1. This observation also suggests that the failure of the iteration in section 2 is likely to imply 
chaos. The identification of particular broken tori provides a starting point for investigating the character and strength of the 
resulting stochasticity. 

The example in section 5 shows that the numerical KAM procedure reproduces many of the same features as the more 
traditional empirical determination using Lyapunov exponents. However, characterising phase space using Lyapunov exponents 
is both time consuming and numerically challenging. They are also difficult to connect to the underlying dynamics for several 
reasons. First, the computation of the Lyapunov exponents requires a numerical solution of the equations of motion. If the orbit 
is chaotic, the initial conserved quantities drift. Secondly, once the orbit enters the vicinity of the homoclinic trajectory that 
has been broadened by nearby resonances, the initially exponentially diverging orbits may appear to converge again. In short, 
orbits with initially positive Lyapunov exponents fall prey to the orbit’s proclivity to explore phase space and locate and get 
stuck near regular islands. So Lyapunov exponents are not ideal for classifying regularity in some arbitrary potential owing to 
the appearance of islands of regularity in the stochastic layers around primary resonances. The same problem affects Poincare 
surface-of-section (SOS) plots; for a single orbit, the SOS plot graphically represents the tendency for weakly chaotic orbits to 
both diffuse to and from regular regions. The numerical KAM approach does not suffer from diffusion and is computationally 
much faster that the standard Lyapunov exponent algorithms for studying the phase space with the same number of initial 
conditions. However, the numerical KAM method requires a rapidly converging Fourier series for accuracy, and this may not 
be satisfied for some problems of astronomical interest. For example, the example in section 5 illustrates that stochasticity 
for very eccentric trajectories in a solar system requires expansion to very high order. On the other hand, the applications 
described in Weinberg (2015a) and Weinberg (2015b)— a bar perturbation to a disc and disc perturbation to a dark-matter 
halo, respectively—will be well described by a small number of terms. 

In addition to analytic estimates of secular evolution and dynamical stability, action variables are widely employed to 
represent phase-space distributions by making use of Jeans’ theorem (Binney & Tremaine 2008). Given a time-independent 
action-space distribution, the positions and velocities of all orbits are known for all time. For many applications, the action- 
angle method will remain a useful and mostly accurate construct. The results of this paper, combined with the galaxy 
applications described in (alias?) and (alias?), suggest that systems with strong deviations from the classical regular systems 
(spherical or triaxial Stackel potentials, two-dimensional axisymmetric disks, etc.) may give rise to dynamics strongly affected 
by, if not dominated by, irregularity. For example, one would not use Jeans’ equations to describe the dynamics in the vicinity 
of Jupiter; (alias?) and (alias?) describe analogous situations for galactic dynamics. 

A secondary goal of this paper is an extension of the well-established tools of secular perturbation theory to the study 
of stochastic dynamics in the context of galaxy evolution. An important motivating application not addressed here or in 
the companion papers is the extension of this technique to the prediction of long-term evolution. The might be done, for 
example, by using the numerical KAM procedure to identify broken tori and then computing diffusion coefficients in the 
initial action space following the numerical approach of Chirikov (1979) or something similar. Phase-space evolution might 
then be followed by a combination of the numerical KAM analyses and a Fokker-Planck type evolution. While feasible, this 
would be a computationally expensive project. 


6.3 Implications for n-body simulations 

Some will no doubt argue that the applicability of the work described here should be established in n-body simulations 
before proceeding further. This gives rise to the question of whether direct n-body simulation suffices to obtain the dynamics 
described in this paper. This is a separate project in itself but in aid of this goal, I offer the following issues for consideration. 

First, it is well established that stochastic layers exist around primary resonances in Hamiltonian systems with multiple 
degrees of freedom. The width of these layers depend on the strength and structure of the perturbation. These layers are 
expected to have lower dimensionality that the phase space itself. Therefore, it is possible that small time-steps, much smaller 
than typically used in n-body simulations, might be necessary for individual particle orbits to resolve (or “feel”) the chaos. 
In other words, large time-steps may be giving the appearance of more regularity than in nature. Anecdotal evidence for this 
issue is found in the computation of Lyapunov exponents. Large time-steps, say 20-100 divisions per characteristic orbital 
time, are too large to capture the effect of exponential growth owing to resonance overlap; time-steps 10-100 times smaller 
were required for convergence of the exponents. 

Secondly, adaptive Poisson solvers generate force fluctuations. Their scale depends on the method. For direct-summation 
codes and tree codes, the spatial scale will be the smoothing length. For PIC-type solvers, the scale will be the cell size. For 
orthogonal-basis codes, the scale is set by the maximum order of the basis functions. The concern is that scattering from these 
fluctuations may interfere with, change or otherwise wash out the underlying Hamiltonian chaos. Some will argue that galaxies 
are naturally noisy, and therefore, if the numerical noise in the force computation washes out the Hamiltonian chaos, then 
it is probably not important. This argument is circular. The existence of astronomical noise does not make n-body particle 


MNRAS 000, 1-25 (2015) 


22 Weinberg 


noise a good stand in. Even in the unlikely case that particle noise emulated the true fluctuations for some application, the 
simulator would need to compute the auto-correlation time for the noise fluctuations to set the correct dynamical time-step. 

The natural noise in dark-matter halos will be a mixture of particle-like and non-particle like. Particle-like noise will come 
from dark substructure and dwarf satellites, molecular clouds and star clusters, to name a few. However, it is likely that these 
are fewer in number than simulation particles and with different distributions than the dark matter or disk particles, and 
therefore the resulting noise spectrum will be different. Non-particle like noise will be generated by the wrapping dark-matter 
streams, continuing gas accretion and other long-term evolutionary processes that result from galaxy assembly. In short, all 
of these processes require precise study before their contribution to the dynamics is well-understood. 


6.4 Topics for future work 

Chaotic dynamics may give rise to anomalous acceleration or scattering caused by overlapping patterns. In the example from 
section 4.1, the stochastic layer that forms around the separatrlx has multiple embedded Islands. Trajectories move through 
the stochastic layer feeling the intermittent effects regular and irregular motion. For example, in computing the Lyapunov 
exponents for section 5 we noticed the tendency for trajectories with positive exponents at early times begin to converge toward 
zero. This is consistent with the behaviour of Hamiltonian intermittency. This well-established behaviour is called stickiness. 
When chaotic trajectories approach the regular regions they stick to their border inducing long periods of almost regular 
motion. This intermittent behaviour determines the magnitude of diffusion. For stronger perturbations, chaotic trajectories 
may appear to suffer from noticeable anomalous kicks and apparently sudden changes in their direction or orientation; such 
effects are related to “radial migration”. 

In cases where the Lyapunov exponents are positive but very small, the trajectories slowly diffuse into the thin chaotic 
layer dehned by the dominant and secondary term(s) in the action-angle expansion of the perturbation. Additional terms 
may yield a complicated network of higher-order resonances, often sticking for very long times to the boundaries of islands 
constituting the so-called edge of chaos regime (Tsallis 2009). This seems less likely to be important to time-limited systems 
such as galaxies. A recent paper by Manos & Machado (2014b) notes that the measured fraction of regular orbits in a barred 
galaxy simulation appears to increase with time. Perhaps this is a consequence of intermittency. 
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APPENDIX A: A COMPUTATIONALLY MORE EFFICIENT ALGORITHM 

The development of the nKAM algorithm in section 2 uses a Fourier series in the new action-angle variables of the perturbed 
Hamiltonian (J,s). While elegant, this set of variables complicated the solution for the new variables in terms of the old 
variables (I, w). Rather than solve equation (11) by a Fourier series in s as described by equation (17), we may write the 
expansion in terms of the old angles w: 

VFi(J,w) = ^t^7i(J)e"*‘'^ (Al) 

1 

where J is the new action vector, w is the old angle vector, and 1 = is a vector of integers whose values are 

1^,1 G 0,lmax for some Imax > 0 as before. This mixed set of variables is more natural for the Type 2 generating function 
(Goldstein, Poole & Safko 2002) but our main motivation here is computational efficiency. This mixed expansion series results 
in iterative algorithm that does not require implicit solution of angle variables. The revised derivations otherwise follows the 
same steps that lead to equation (9). 

To see this new formulation facilitates the implicit solution of the new canonical variables as a function of the old, write 
the expansion of the right-hand side of equation ( 6 ) to first order in the perturbation Wi — 0[Hi] explicitly: 

/ dwC^^'' \ f 91Ti 1 9 

=IiIj+ = I Ho(J) + fl(J) • -h I + C>(W ) (A2) 

where we have used the notation f2j(J) = dHo/dljjj for compactness. If we subtract the terms in {} from both sides of 
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equation (A2), after rearranging and grouping the new terms, we have: 

dWi 


H„e»(J) - ii'o(J) - n(J) 


dw 




Ho J + 


Hi J + 


dWi 

9w 

dWi 

dw 


- Ho(J) - n(J) 


-Hi(J,w) 


dWi 


dw 


+ 


(A3) 


The left-hand side is the difference between Hnsw{J) and the first-order Taylor series. Terms on the right-hand side of 
equation (A3) are second order in the perturbation Vi. As in the section 2.2.1, we can exploit this to determine the first-order 
contribution once and for all, and iterate to improve the second-order contribution. 

We proceed as follows. Consider the action-angle expansion of equation (A3) in the Fourier series defined by equation 
(Al). The first two terms on the left-hand side are constants with respect to w and therefore only contribute to coefficients 
ti7i with 1 = 0. From equation (3), these terms make no contribution to the new action, and from equation (4), these terms 
only result in a constant offset in new angle. As in section 2, we see from equations (3) and (4) that the 1 = 0 transform of 
equation (A3) generates an uninteresting phase shift in the solution. The Fourier coefficient of the third term on the left-hand 
side of equation (A3) is —il ■ fl(J)'n7i(J) while that for the fourth term is determined once and for all by direct quadrature. 
Since the right-hand side is second order in Wi , we can immediately solve for the first-order contribution to zui by demanding 
that the first-order contribution to the LHS of equation (A3) vanishes: 


411 


(J) = 


dwHi(I, w)e 


(A4) 


1- n(J) (27r)3 

where I = J and s = w to first order. Again, we can safely ignore the 1 = 0 which yields an constant offset in s by equation 
(4). We may now substitute this first-order solution 


W™(J,w) = ^t:u,(J)We*''"' 


(AS) 


into equation (A3) to get: 
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(J) + 


1- r2(J) (27r)3 
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w j — Hi(J, w) 


(A6) 


Equation (A6) may be similarly iterated to obtain 
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dw 
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dw 


dw 
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(A7) 


or using H — Hq + Hi, 
^1"’ (J) 


(J) -I- 


i 1 

1- n(J) (27r)3 


^ dw 


H J + 


dWf- 

dw 


- H(J,w) - fl 


dWi 

dw 


[n-4-| 


(AS) 

The iteration may be terminated when the value of maxi Iwj”' — is sufficiently small. 

Equations (A4), (A6), (A7), and (AS) require angle integrals over a range of 1. The restricted three-body problem 
example in section 5 and the galactic-bar perturbation in (alias?) considers a two-dimensional orbital distribution with a 
non-axisymmetric perturbation. For the application in (alias?), we will consider azimuthally symmetric perturbations only; 
Lz will remain constant for all orbits and only terms zui with I 3 = 0 will be non-zero. Therefore, our angle integrals in the 
equations above are two dimensional for this case as well. By restricting our series to |Zi| G [0, li,max] and IZ 2 I G [0, l 2 ,max], the 
angle integrals in equations (A4), (A6) and (A7) have the same form as a discrete Fourier transform and may be efficiently 
performed using the FFT algorithm. The two-dimensional restriction is for numerical convenience only; see (alias?) for a 
discussion of the full three-dimensional perturbation. 


APPENDIX B: GRID CHOICE AND INTERPOLATION 

As described in section 2.4.1, we divided our action grids into three classes: 1) two-dimensional problems with action space 
(7 i,72 ); 2) three-dimensional problems with action space (IijHjIa) with one action conserved, I3 = Lz say; and 3) the fully 
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general three-dimensional problem. Our software implementation includes all three classes but only (1) and (2) are used in 
this and the companion papers. In general, the distribution of I will not be on a regular grid, and therefore interpolation and 
evaluation of partial derivatives require specialised methods. We have explored spline fits, interpolating simplices (Riachy, 
Mboup & Richard 2011) and the Shepard method (Shepard 1968). Numerical experiments suggests that the Shepard method 
is the most robust and generally accurate on test problems. The authors of ALGLIB (www.alglib.net) suggest that the radial 
basis function (RBF) methods (e.g. Buhmann 2003) further improve the performance of Shepard’s algorithm, but this has 
not been investigated here. 

Our implementation of Shepard’s algorithm begins with an irregular point set denoted as li for i = 1,... ,n. Let the 
desired field F(I) evaluated on this set be denoted as Fi. The original Shepard approximating function (Shepard 1968) is 
defined as 


F(I)=^ru.(I)F. 

i=l 


where the weight functions Wi are 


«;i(I) = 


H[5-di{I)Y 

j:UH[s-dAi)Y 


with 5 > 0, di = ||I — Ii||, the standard Euclidean norm, and H[-] defined by 


m 


q, q^O 
0, q < 0 


(Bl) 


(B2) 


(B3) 


The smoothing parameter 5 and the exponent r are tuning constants. Generally, if exponent r has a value in [0,1] then F 
has peaks at the grid points, and if r > 0, the approximating function tends toward a well defined slope at the nodes. In this 
paper, we use the modified Shepard method implementations by Renka (1988a,b). This updated algorithm generalises the 
simple weights of the original algorithm by augmenting the approximation of Fi by a quadratic multinomial whose coefficients 
are determined by the weighted least-squares fit to the N nearest neighbours to I. The weights are 

H[n{NA - di{V)Y 


«.(I) = ■ 


[niNAdYl)? 


(B4) 


where ri{Nq) is the radius about li containing Nq points. We used the recommended default value of Nq = 13 and 17 for 
the two- and three-dimensional cases from Renka (1988a,b). The derivatives (e.g. for gradients of field quantities) are easily 
computed analytically from the quadratically augmented form of Fi. Finally, at each node in the grid, we store the full orbit 
computation, tabulated as function of angle. These are linearly as needed to produce a full phase-space evaluation at (I,w). 
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